  function drawsmat=genmultrand(distvec,avec,bvec,nd,lowerBound,upperBound) 
% function drawsmat=genmultrand(distvec,avec,bvec,nd); 
%
% GENerate MULTiple RANDom draws from a density dist with 
% location coefficients a, b. 
% For a single set of draws use genrand.m 
% 
% Inputs 
% ========
% dist     G,N,B,W,I,U 
% a and b  For G and B, alpha and beta  
%          For N density mean and sd (NOT the variance!);
%          For I and W, v( degrees of freedom) and s (scale)
%          For U, limits [a,b] 
%          For L, create a linspace of nd points equally spaced between 
%              [a,b]
%
% nd       (optional) if not provided nd = 20,000 
%
% Created               6/10/04
% Last Modified         1/22/05  Added L as an option 
% Alejandro Justiniano  ajustiniano@imf.org 
% ------------------------------------------------------------
if nargin < 4 
    nd=20000; 
else
    if nd < 0 
        error('Need positive number of draws')
    end
end   

avec=avec(:); 
bvec=bvec(:); 
np=size(avec,1); 
drawsmat=zeros(nd,np); 

for jj=1:np
    jj;
    np;
    dist=upper(char(distvec(jj)));
    a=avec(jj); 
    b=bvec(jj); 
    lb=lowerBound(jj); 
    ub=upperBound(jj); 
    
    if isstr(a) ==1 | isstr(b) == 1
        error('Inputs a or b must be scalars') 
    end 
    if ( ( strcmp(dist,'N')==0 )  && ( strcmp(dist,'U')==0 ) && ( strcmp(dist,'C')==0 ) )  ; 
        if a < eps/100; error('Location 1 cannot be negative'); end;  
        if b < eps/100; error('Location 2 cannot be negative'); end;   
    end
    switch dist 
        case 'N'
            if b < (eps/100); 
                error('Variance of a normal cannot be negative'); 
            end;  
            draws=a+b*randn(nd,1);  
        case 'G' 
            draws=gamrnd(a,b,nd,1); 
        case 'B' 
            draws=betarnd(a,b,nd,1);
        case 'U' 
            draws=unifrnd(lb,ub,nd,1);
        case 'W' 
            draws=gamrnd(a/2,2/b,nd,1);
            draws=1./draws; 
        case 'I' 
            draws=gamrnd(a/2,2/b,nd,1);
            draws=1./draws; 
            draws=sqrt(draws);
        case 'L' 
            draws=linspace(a,b,nd); 
        case 'C' 
            if b ~=0; 
                error('Calibrated parameter must have zero variance'); 
            end; 
            disp('Calibrated parameter....') 
            draws=a*ones(nd,1); 
    end
    drawsmat(:,jj)=draws; 
end; 